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We study the effects of magnetic fields on the evolution of differentially rotating neutron stars, 
which can be formed in stellar core collapse or binary neutron star coalescence. Magnetic braking 
and the magnetorotational instability (MRI) both act on differentially rotating stars to redistribute 
angular momentum. Simulations of these stars are carried out in axisymmetry using our recently 
developed codes which integrate the coupled Einstein-Maxwell-MHD equations. We consider stars 
with two different equations of state (EOS), a gamma-law EOS with V = 2, and a more realistic 
hybrid EOS, and we evolve them adiabatically. Our simulations show that the fate of the star 
depends on its mass and spin. For initial data, we consider three categories of differentially rotat- 
ing, equilibrium configurations, which we label normal, hypermassive and ultraspinning. Normal 
configurations have rest masses below the maximum achievable with uniform rotation, and angular 
momentum below the maximum for uniform rotation at the same rest mass. Hypermassive stars 
have rest masses exceeding the mass limit for uniform rotation. Ultraspinning stars are not hyper- 
massive, but have angular momentum exceeding the maximum for uniform rotation at the same rest 
mass. We show that a normal star will evolve to a uniformly rotating equilibrium configuration. An 
ultraspinning star evolves to an equilibrium state consisting of a nearly uniformly rotating central 
core, surrounded by a differentially rotating torus with constant angular velocity along magnetic 
field lines, so that differential rotation ceases to wind the magnetic field. In addition, the final state 
is stable against the MRI, although it has differential rotation. For a hypermassive neutron star, the 
MHD-driven angular momentum transport leads to catastrophic collapse of the core. The resulting 
rotating black hole is surrounded by a hot, massive, magnetized torus undergoing quasistationary 
accretion, and a magnetic field collimated along the spin axis — a promising candidate for the central 
engine of a short gamma-ray burst. 

PACS numbers: 04.25.Dm, 04.30.-w, 04.40.Dg 



I. INTRODUCTION 

Differentially rotating neutron stars can form from the 
collapse of massive stellar cores, which likely acquire 
rapid differential rotation during collapse even if they 
are spinning uniformly at the outset 0,0 (see also 0)- 
Differential rotation can also arise from the mergers of 
binary neutron stars 0, |H El- In these new-born, dy- 
namically stable, neutron stars, magnetic fields and/or 
viscosity will transport angular momentum and cause 
a substantial change in the configurations on a secular 
timescale. 

Some newly-formed differentially rotating neutron 
stars may be hypermassive. Specifically, the mass lim- 
its for non-rotating stars [the Oppenheimer-Volkoff limit] 
and for rigidly rotating stars (the supramassive limit, 
which is only about 20% larger) can be significantly 
exceeded by the presence of differential rotation U- 
Mergers of binary neutron stars could lead to the for- 
mation of such hypermassive neutron stars (HMNSs) as 
remnants. This possibility was foreshadowed in New- 
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tonian post-Newtonian and in full general rela- 
tivistic simulations 0]. The latest binary neutron star 
merger simulations in full general relativity 0, 0, 0] 
have confirmed that HMNS formation is indeed a possi- 
ble outcome. HMNSs could also result from core collapse 
of massive stars. 

Differentially rotating stars tend to approach rigid ro- 
tation when acted upon by processes which transport an- 
gular momentum. HMNSs, however, cannot settle down 
to rigidly rotating neutron stars since their masses exceed 
the maximum allowed by rigid rotation. Thus, delayed 
collapse to a black hole and, possibly, mass loss may re- 
sult after sufficient transport of angular momentum from 
the inner to the outer regions. Several processes can act 
to transport angular momentum and drive the HMNS 
to collapse. Previous calculations in full general relativ- 
ity have modeled the evolution of HMNS driven by vis- 
cous angular momentum transport |12| and by ang ular 
momentum loss due to gravitational radiation [1(1 ]. In 
both cases, the core of the HMNS eventually collapses to 
a black hole. However, in the case of viscosity-driven 
evolution, a large accretion torus is found to develop 
around the newly-formed black hole, while for gravita- 
tional wave-driven evolution, the disk present after col- 
lapse is very small. The size of the disk is of crucial im- 
portance, because if a large disk is produced, the post- 
collapse system may produce a short-duration gamma- 
ray burst (GRB). 
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The merger of binary neutron stars has been pro- 
posed for many years as an explanation of short-hard 
GRBs 0,13- Possible associations between short GRBs 
and elliptical galaxies reported recently [l5j make it un- 
likely that short GRBs are related to supernova stellar 
core collapse. The merger of compact-object binaries 
(neutron star-neutron star or black hole-neutron star) is 
now the favored hypothesis for explaining short GRBs. 
According to this scenario, after the merger, a stellar- 
mass black hole is formed, surrounded by hot accretion 
torus containing ~ 1-10% of the total mass of the system. 
Energy extracted from this system, either by MHD pro- 
cesses or neutrino-radiation, powers the fireball for the 
GRB. The viability of this model depends on the pres- 
ence of a significantly massive accretion disk after the 
collapse of the remnant core, which in turn depends on 
the mechanism driving the collapse. 

Though magnetic fields likely play a significant role 
in the evolution of HMNSs, the numerical tools needed 
to study this problem have not been available until re- 
cently. In particular, the evolution of magnetized HMNSs 
can only be determined by solving the coupled Einstein- 
Maxwell-MHD equations self-consistently in full general 
relativity. Recently, Duez et al. |l6j and Shibata and 
Sekiguchi yjj independently developed codes designed 
to do such calculations for the first time (see also [HI). 
The first simulations of magnetized hypermassive neu- 
tron star collapse (assuming both axial and equatorial 
symmetry) were reported in |l9j |. and the implications 
of these results for short GRBs were presented in [20| . 
These simulations proved that the amplification of small 
seed magnetic fields by a combination of magnetic wind- 
ing and the magnetorotational instability (MRI) is suf- 
ficient to trigger collapse in hypermassive stars on the 
Alfven timescale, confirming earlier predictions 0, l2l| . 
In the present work, we describe these collapse calcula- 
tions in more detail. 

We also compare the results for hypermassive stars 
with the evolution of two differentially rotating models 
below the supramassive limit in order to highlight the 
qualitatively different physical effects which arise in the 
evolution. Given a fixed equation of state (EOS), the 
sequence of uniformly rotating stars with a given rest 
mass has a maximum angular momentum J max - A non- 
hypermassive star having angular momentum J > J max 
is referred to as an "ultraspinning" star. We perform sim- 
ulations on the MHD evolution of two nonhypermassive 
stars - one is ultraspinning and the other is not; we refer 
to the later as "normal." Instead of collapsing, they set- 
tle down to a new equilibrium state after several Alfven 
times. The normal star settles down to a uniformly ro- 
tating configuration. In contrast, the ultraspinning star 
settles down to a nearly uniformly rotating central core, 
surrounded by a differentially rotating torus. In this new 
equilibrium, the system has adjusted to a state where the 
angular velocity is constant along the magnetic field lines, 
which means that the residual differential rotation ceases 
further magnetic winding. In addition, we find that the 



final state is also stable against the MRI, although it has 
differential rotation. 

The key subtlety in all of these simulations is that the 
wavelengths of the MRI modes must be well resolved on 
the computational grid. Since this wavelength is pro- 
portional to the magnetic field strength, it becomes very 
difficult to resolve for small seed fields. However, the 
simulations reported here succeed in resolving the MRI. 

The structure of this paper is as follows. We give an 
overview of the MHD effects acting on differentially ro- 
tating stars in Sec. [H] The initial models are briefly 
discussed in Sec. IIIII In Sec. IIVI we summarize the set 
of coupled Einstein-Maxwell-MHD equations which are 
solved during the simulations. We outline our numeri- 
cal methods for the evolution in Sec. EJ and present our 
simulation results in Sec. lVIl Finally, we summarize and 
discuss our main results in Sec. I VI II In what follows, we 
assume geometrized units such that G = c = I. 



II. OVERVIEW OF MHD EFFECTS 

Two distinct processes which are known to trans- 
port angular momentum in differentially rotating mag- 
netized fluids are magnetic braking 0, El H3, and 
the MRI 0, H3 • Magnetic braking tr ansp orts angular 
momentum on the Alfven time scale 0, |21| : 
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where R is the radius of the HMNS and va is the Alfven 
speed. 

At early times, the effects of magnetic braking grow 
linearly with time. This can be seen by considering the 
magnetic induction equation in a perfectly conducting 
(MHD) plasma [see Eq. (gSJl below]: 



where 



dtB 1 + d J (v°B i -v l B j ) = 



B* = V7£T = ^n v F* 



(2) 



(3) 



In the above formulae, 7 is the determinant of the spatial 
metric, n M = (—a, 0, 0, 0) is the normal to the spatial hy- 
persurface, a is the lapse, F*^ v is the dual of the Faraday 
tensor, and v l = u l /u° is the 3-velocity of the fluid. Now, 
if the magnetic field is weak and has a negligible back- 
reaction on the fluid, the velocities will remain constant 
with time. In cylindrical coordinates, we have (assuming 
axisymmctry) 

d t B l « , (i = w,z) (4) 

d t B v = -d^B* -v v B l ) 

« d^B 1 ) , {i = m,z) (5) 

where w is the cylindrical radius, and where we have 
used the fact that w ro = v z — at t = and remains 
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so under these assumptions. Then, since v v = CI (the 
angular velocity), 

d t B v re B^n + fld.B' (i = zu,z). (6) 

The second term vanishes by Maxwell's equations (the 
no-monopole constraint, diB 1 — 0) and the assumption 
of axisymmetry (d v — 0). At early times, Eq. © indi- 
cates that the toroidal component of the field B T = wB^ 
grows linearly according to 

B 7 (t; w, z) ps twB l (0; w, z)difl(0; zu, z) [i = zu, z) . 

(7) 

The growth of B T is expected to deviate from this lin- 
ear relation when the tension due to the winding up of 
magnetic field lines begins to change the angular velocity 
profile of the fluid. 

The MRI is present in a weakly magnetized, rotating 
fluid wherever d^fl < [H El When the instability 
reaches the nonlinear regime, the distortions in the mag- 
netic field lines and velocity field lead to turbulence. To 
estimate the growth timescale £mri and the wavelength 
of the fastest growing mode A max , we make use of a sim- 
ple Newtonian linear analysis given in j2|J (see also [27j ) ■ 
Linearizing the MHD equations for a local patch of a ro- 
tating fluid and imposing e I ' k ' x "" i ' dependence on the 
perturbations leads to the dispersion relation given in 
Eq. (125) of [2g. Specializing this equation for a constant 
entropy star and considering only modes in the vertical 
direction, this reduces to 

w 4 -[2(k-v A ) 2 + K 2 ]w 2 + (k-v A ) 2 [(k- Vj 4) 2 + K 2 -4ri 2 ] = , 

(8) 

where = B / ^Airp is the (Newtonian) Alfven velocity, 
p is the mass density and k is the "epicyclic frequency" 
of Newtonian theory: 
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(9) 



We consider vertical modes (k = ke z ) since we are only 
looking for an estimate of A max , and since these are likely 
to be the dominant modes. 

Since the dispersion relation in Eq. JHJ) is quadratic in 
lo 2 , it can be easily solved for to 2 and then minimized to 



find the frequency of the fastest-growing mode, ui n 

s 1 

4(s 2 + k 2 ) ~ 4 



9 In r 



(10) 



where s 2 = 4f2 2 — n 2 
growth rate corresponds to 



-2Qdfl/dhir. This maximum 
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For the growth time and wavelength of the fastest grow- 
ing mode, we then have 

-l 



£mri 



2{dfl/d\nzu 
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? + 2k 2 



(12) 
(13) 



In order of magnitude, 
Amax ~ 2irv z A /ri 
~ 3 cm 



4000 rads"iy \ 10 12 G 
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4000 rad s~l 



(14) 



(15) 



Here ft = 4000 rad s" 1 corresponds to a rotation pe- 
riod P = 1.57 ms. For realistic HMNS magnetic fields, 
Amax will be much smaller than R. We note that, since 
Amax <x va, larger magnetic fields will result in longer 
MRI wavelengths. When A max ^ R, where R is the 
equatorial radius of the star, the MRI will be suppressed 
since the unstable perturbations will no longer fit inside 
the star. This is why the MRI is regarded as a weak-field 
instability. Typically, we set magnetic field amplitudes 
so that A max ~ i?/10 for the models we consider here. 
We note that turn is independent of the strength of the 
seed magnetic field. The MRI always grows on a dy- 
namical timescale for a sufficiently differentially rotating 
configuration. Hence, the MRI is likely to be very im- 
portant during the early evolution for realistic HMNSs. 
However, the resulting angular momentum transport is 
governed by the turbulence and is thus expected to occur 
on a timescale longer than iMRi- 

Magnetic fields and turbulence tend to transport spe- 
cific angular momentum from the rapidly rotating inner 
region of a differentially rotating star to the more slowly 
rotating outer layers. This causes the inner part to con- 
tract and the outer layers to expand. Since hypermassive 
stars depend on their strong differential rotation for sta- 
bility, this angular momentum transport process likely 
leads to collapse. However, in Section IVI CI we show 
that very different behavior can result for rapidly rotat- 
ing nonhypermassive models. In the example we explore, 
the star readjusts to a new equilibrium state consisting 
of a nearly rigidly rotating core surrounded by a differen- 
tially rotating torus in which the magnetic field lines are 
everywhere orthogonal to the gradient of the angular ve- 
locity (i.e., B^djfl = 0). Hence, magnetic winding shuts 
down even though the configuration is still differentially 
rotating. This possibility has been discussed previously 
by Spruit |2^| in the context of Newtonian theory. 



III. INITIAL MODELS 

In order to study the effects of rotation and EOS, 
we evolve four representative differentially rotating stars, 
which we call "A" , "Bl" , "B2" and "C" . Their properties 
are listed in Table [I] Stars A and C are hypermassive; 
stars Bl and B2 are not. These configurations are all 
dynamically stable. 

Stars A, Bl and B2 are constructed using a T — 2 poly- 
tropic EOS, P — KpQ, where P, K, and po are the pres- 
sure, polytropic constant, and rest-mass density, respec- 
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TABLE I: Initial Models 



Case 


EOS 


Mo/M ,TOv a 


M /M OiSup b 


M/M snp ° 


i? oq /M d J/M 2 c 


T Iot /\W\ l 




P c /M h 


A 


r = 2 


1.69 


1.46 


1.49 


4.48 


1.0 


0.249 


0.33 


38.4 


Bl 


r = 2 


0.99 


0.86 


0.89 


8.12 


1.0 


0.181 


0.40 


103 


B2 


r = 2 


0.98 


0.85 


0.86 


4.84 


0.38 


0.040 


0.34 


105 


C 


hybrid 


1.28 


1.14 


1.17 


2.75 


0.82 


0.241 


0.185 


15.5 



a The ratio of the rest mass Mo to the TOV rest mass limit for the given EOS. 

b The ratio of the rest mass Mo to the rest mass limit for uniformly rotating stars of the given EOS (the 
supramassive limit). If this ratio is greater than unity, the star is hypermassive. 

c The ratio of the ADM mass M to the gravitational mass limit for uniformly rotating stars of the given EOS. 

d The equatorial coordinate radius R cq normalized by the ADM mass. 

e The ratio of the angular momentum J to M 2 (the angular momentum parameter) . 

f The ratio of the rotational kinetic energy to the gravitational binding energy [see Eqs. 15611 and 1581 1. 

s The ratio of the angular velocity at the equator to the central angular velocity. 

h The initial central rotation period P c normalized by the ADM mass. 



tively. (In [Tj] , which considered evolution with shear vis- 
cosity, star A was referred to as "star I" and star Bl was 
referred to as "star V"). The rest mass of star A exceeds 
the supramassive limit by 46%, while the rest masses of 
stars Bl and B2 are below the supramassive limit. The 
angular momentum of star Bl exceeds the maximum an- 
gular momentum (J max ) for a rigidly rotating star with 
the same rest mass and EOS, whereas star B2 has angular 
momentum J < J max . Thus, star Bl is "ultraspinning," 
while star B2 is "normal." Stars A, Bl and B2 may 
be scaled to any desired physical mass by adjusting the 
value of K |2^. In general, M oc K n / 2 , where n is the 
polytropic index (r = 1 + 1/n, here n = 1). For example, 
choosing K = 2.42 x 10 5 g -1 cm 5 s -2 gives the maximum 
ADM mass (rest mass) of 2.12M© (2.32M Q ) for spherical 
neutron stars and 2.42M Q (2.67M©) for rigidly rotating 
neutron stars. 

In order to consider the effects of a more realistic neu- 
tron star equation of state, star C is constructed from a 
cold hybrid EOS defined as follows: 



P = P, 



cold 



Kxfc for PO < Pr, 
K 2 p L 2 for PO > Pu 



(16) 



We set r x = 1.3, T 2 = 2.75, K x = 5.16 x 10 14 cgs, 
K 2 = K x( ^- V \ and p nuc = 1.8 x 10 14 g/cm 3 . 
With this EOS, the maximum ADM mass (rest mass) 
is 2.01Af Q (2.32M©) for spherical neutron stars and 
2.27M Q (2.6OM0) for rigidly rotating neutron stars, 
which are similar values to those in realistic stiff 
EOSs 29]. Star C exceeds the supramassive limit by 
14%. The various parameters of star C are chosen in or- 
der to more closely mimic the HMNSs formed through 
binary neutron star mergers with realistic equations of 
state in [To|. 

Following previous papers (e.g, 0, 0, IH, HU), we 
choose the initial rotation law u°u v = A 2 (fl c — ft), where 
is the four-velocity, f2 c is the angular velocity along 



the rotational axis, and = u v /u° is the angular veloc- 
ity. In the Newtonian limit, this rotation law becomes 



a. 



(17) 



The constant A has units of length and determines the 
steepness of the differential rotation. In this paper, A 
is set equal to the coordinate equatorial radius R eq for 
stars A, Bl, and B2, while A = 0.8i? cq for star C. The 
corresponding values of Q cq /Q c are shown in Table [fl 
(where Q cq is the angular velocity at the equatorial sur- 
face). The magnitude of the angular momentum is seen 
from the Kerr parameter a/M = q = J/M 2 . Stars A, Bl, 
and C have q = 1.0, 1.0, and 0.82, respectively. These 
stars rotate very rapidly and are highly flattened due to 
centrifugal force. Star B2, on the other hand, has a com- 
paratively low angular momentum parameter: q = 0.38. 

We must also specify initial conditions for the magnetic 
field. We choose to add a weak poloidal magnetic field 
to the equilibrium model by introducing a vector poten- 
tial of the following form A v = va 2 max.[Ai,(P — P cu t), 0], 
where the cutoff P C ut is 4% of the maximum pressure, and 
Ai, is a constant which determines the initial strength of 
the magnetic field. We characterize the strength of the 
initial magnetic field by C = max(6 2 /P), i.e. the maxi- 
mum value on the grid of the ratio of the magnetic en- 
ergy density to the pressure. We choose At such that 
C ~ 10 _3 -10 -2 . We have verified that such small ini- 
tial magnetic fields introduce negligible violations of the 
Hamiltonian and momentum constraints in the initial 
data. 
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IV. BASIC EQUATIONS 

A. Evolution of the gravitational fields 

We evolve the 3-metric 7^ and extrinsic curvature 
Kij using the Baumgarte-Shapiro-Shibata-Nakamura 
(BSSN) formulation The fundamental variables for 
BSSN evolution are 



= 


— ln[det(7<j)] , 


(18) 


lij = 


e^Hi , 


(19) 


K = 


l lJ K l3 , 


(20) 


Aij = 




(21) 


r = 


-f\ 3 (or Pi = ^'*7y, fc ) . 


(22) 



The evolution equations for these variables are as follows: 



(dt - Cp)jij = -2aAij 
(d t -£p)<t> = ~aK 

(d t -C p )K -- 



7^Dj DiO. + -aK 2 



+aA lj A lj + Ana(p + S) 

DiDja + ( 
-a(KAij -2A a A l j ) 



e- 4 * (-DiDja + a(R i3 - 8ttS %3 )) 



(23) 
(24) 

(25) 

TF 

(26) 



where a is the lapse function, (3 l is the shift, Cp is the 
Lie derivative along the shift, and Di is the covariant 
derivative with respect to the spatial 3-metric. The Ricci 



tensor R i3 can be written as the sum 



Rij — Rij + R- . 



(27) 



Here Rf- is 



-2A-Dj0 - 2%D l D^ 

+4(D i <f>)(D j <j>) - 4 7y -(^0)(A0), (28) 



where D l = j^Dj. The "tilde" Ricci tensor R i3 is the 
Ricci tensor associated with 7^ , and is computed by 



R, 



\l lm lHM + lk(idj)f k + f k t {ij)k + 



2f?/,r 



l(i L j)km + F im Tklj 



where 



Fijk = 7>{lH,k +lik,j - ljk,i) ■ 



(29) 



(30) 



The evolution of f % is given by 
d t f l = d j (2aA ij +Cpj ij ) 



1 



l 0k p l ,3k + fl l ^\k 3 -VP\ 3 

+^r l p\ 3 + ^f\ J 



(31) 



2A %3 dja 



-2a [ - QA l3 <j) :3 - r jk A jh + 871-7" 5,- 



The matter source terms p, Si 3 , Si and S are related to 
the stress-energy tensor T^ v as follows: 



P = 

S t = 

Sij — 

s = 



n a npT af3 , 
S\ . 



(32) 



In the code of Duez et al. [l(|. additional constraint 
damping terms are included in the BSSN evolution sys- 
tem as described in jH Hi (see Eqs. (45) and (46) 
of Eqs. (6)-(8), (11), (13) and (15) of H). The 
gauge conditions used with this code are the following 
hyperbolic driver conditions 03 : 



d t a 
d t A 
d 2 ? 



aA , 

—a\(adtK + a 2 dta 
bi{ad t r -bidtft) , 



a 3 e 



b aK) , (33) 
(34) 



where a±, a 2l 03, 61, and b 2 are freely specifiable con- 
stants. We usually choose a\ = 0.75, b\ = 0.15, a 2 and 
6 2 between 0.34/M and 0.56/M , and a 3 between 0.17/M 
and 0.28/M, where M is the ADM mass of the star. 
For runs with excision, we use a\ = bi = 0.75, 03 = 0, 
0.34/M <a 2 < 0.56/M and b 2 = a 2 . 

In the code of Shibata and Sekiguchi the following 
dynamical gauge conditions are used: 



dta — —aK, 
dtf? = ^(F 3 +AtdtF 3 ), 



(35) 
(36) 



where At is the timestep, and Fj is the function defined 
in Eq. I|22|) . For the evolution, constraint damping terms 
are added to the equations for <j> and Ay to suppress 
high-frequency noise and maintain the accuracy of the 
Hamiltonian constraint and tv(Ai 3 ) = 0. 

B. Evolution of the electromagnetic fields 

The evolution equation for the magnetic field in a 
perfectly conducting MHD fluid (F^ v u v = 0) can be 
obtained in conservative form by taking the dual of 
Maxwell's equation -F] m „.a] = 0. One finds 



1 



-9 



(37) 
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where y/—g = ce^/j, F a ^ is the Faraday tensor, and F* a @ 
is its dual. Using the fact that the magnetic field as mea- 
sured by a normal observer n a is given by B % = n^F*^ 1 , 
the time component of Eq. (|37|l gives the no-monopole 
constraint djB 3 = 0, where B 3 = B 3 . The spatial 
components of Eq. IpFTjl give the magnetic induction equa- 
tion, which can be written as 



d t B l + d 3 {v 3 B i -v l B 3 ) = 



C. Evolution of the hydrodynamics fields 



(38) 



E2 



The evolution equations for the fluid are as follows 0, 



d t p* + dj{p*v 3 ) = , 



(39) 



dtSi + d^ctyfyT'i) = ^a^T^g^i , (40) 

(41) 



dtf + 0i(aVr T = s 



where the density variable is p* = a^/jpou , the 
momentum-density variable is Si = ay^T?, the energy- 
density variable as adopted by Duez et al. fl6l ] is f = 
a 2 ypy T 00 — p» , and the source term s is 

s = -<*/yT' u 'V„n /I 

= aV7 [{T 00 f3 i f3 j + 2T Ql (3 3 + T l3 )K l} 

-(T oa p l + T m )d t a] . (42) 

The MHD stress-energy tensor is given by 

= { Po h + b 2 )u^u v + (P + b 2 /2)g^ - b»b v , (43) 



where 6 M = u v F* vtl /\/4tt, the specific enthalpy is given 
by h = 1 + e + Pj po, e is the specific internal energy, and 
b 2 = b^tt 1 . 

In the code of Shibata and Sekiguchi [l?J, the energy 
evolution variable is chosen to be ^f^n^rivT^ 1 ' = r + p*, 
and the evolution equation may be obtained by adding 
Eq. to Eq. gTJ}. 

The MHD system of equations is completed by a choice 
of EOS for the evolution. For stars A, B2, and B2, we 
adopt a T-law EOS P = (T - l)p e, with T = 2. For 
star C, we adopt the following hybrid EOS: 

P = P CO M + (Ta, - l)A)(e - Scow)- (44) 

Here, Pr.o iH and £ C oid denote the cold component of P 
and e |l7| . The conversion efficiency of kinetic energy to 
thermal energy at shocks is determined by Tth , which we 
set to 1.3 to conservatively account for shock heating. 



D. Diagnostics 

We monitor several global conserved quantities to 
check the accuracy of our simulations. The ADM mass 



M and angular momentum J are defined as integrals over 
surfaces at infinity as follows [35| : 



M = 



Ji 



lOTT J r=oa 



/ r—oo 
k 



8tt 



x j KFd 2 S< 



)d 2 S l ,(45) 
(46) 



In cases for which no singularity is present on the grid, 
these surface integrals can be converted to volume inte- 
grals using Gauss's theorem (see Appendix A of 32]): 



M 



^(Po + - ^-K 2 ) 



■ptjfc-r ... 

167T ^ 



1 



(47) 



16tt 



-R 



<f;< 



T — F- k 



12-K l07T 



(48) 



These integrals should be exactly conserved. However, 
using finite grids, we are unable to perform this integral 
out to infinity, and we expect to see mass and angular mo- 
mentum losses due to outflows (of fluid, electromagnetic 
fields, and/or gravitational waves) through the bound- 
aries. These fluxes can be measured, however, and are 
found to be quite small. 

In axisymmetry, the volume integral for the angular 
momentum (which is entirely in the z-direction) simpli- 
fies considerably |36| : 



J = I S v d 6 x 



(49) 



An additional conserved quantity is the total rest mass 
M : 



M = / Pi d 6 x 
Jv 



(50) 



In axisymmetry, gravitational radiation carries no angu- 
lar momentum, and in this case our GRMHD codes are 
finite differenced such that Mo and J are identically con- 
served in the absence of flux through the boundaries. 
Hence, M and J are not useful diagnostics when volume 
integrals l)49|) and (|50|l are applicable. 

For runs with black hole excision, a volume integral 
must be replaced with an integral over an inner surface 
surrounding the black hole plus a volume integral ex- 
tending over the rest of the grid (see |32( for details). 
The integral for J is then no longer identically conserved 
by our numerical scheme, and the total angular momen- 
tum is only constant to the extent that the excision evo- 
lution is accurate. During excision evolutions, we sep- 
arately track the rest mass and angular momentum of 
matter outside the hole by carrying out the integrals in 
Eqs. H49|) and l|50|l over the region outside the apparent 
horizon. Though no longer exact, these integrals allow 
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us to estimate the rest mass and angular momentum of 
the accretion torus. 

When a black hole is present , w e detect it by using 
an apparent horizon finder (see |3jj for details). As the 
system approaches stationarity, the apparent horizon will 
approach the event horizon. From the surface area of 
the apparent horizon .4 ah, we compute the approximate 
irreducible mass M- 1IT by 



V-4ah/167t 2 



(51) 



In order to check the accuracy of our simulations, we 
monitor the L2 norms of the violation in the constraint 
equations. In terms of the BSSN variables, the constraint 
equations become, respectively, 



= H = -'•'/),/;,." - — R 



(52) 



,50 



= M i = D J (e 64 'A jt ) - -e^D l K - Sne^S 1 .(53) 

O 

We normalize Tt and M. % and compute the L2 norms on 
the grid as described in [38|. 

In order to understand the evolution of the magnetic 
field, it is useful to compute field lines. Below, we plot 
field lines corresponding to the poloidal magnetic field. 
In axisymmctry, these field lines correspond to the level 
surfaces of A v (see Appendix^) , which is computed from 
_B ro and B z . To visualize the toroidal field, we also plot 
the 3D field lines projected onto the equatorial plane (see 
Appendix for details of the method). 

We measure several invariant energy integral diagnos- 
tics during the evolution. We define the adiabatic inter- 
nal energy i?i n t,ad, the internal energy from heat, -Eheat, 
rotational kinetic energy T rot , the electromagnetic energy 
-Eem, and gravitational potential energy W, as follows: 



-Eint,ad 
-Ehcat 
Trot 

W 



(/Ooehcat)dV , 



1 



HT^dV/u , 



v 2 

n^n v T^dV / \au°) 



v 

M - Mn - E, 



^int,ad ^Ticat 



Ei 



EM 



(54) 

(55) 

(56) 

(57) 
(58) 



where dV = au -^fd x is the proper 3- volume element, 
^fluid = Pohu^u" + PgV is the perfect fluid stress-energy 
tensor, = b 2 u^u v + b 2 g» v /2 - b^b y is the stress- 

energy tensor associated with the electromagnetic field, 
£coid refers to a the cold initial polytrope or hybrid EOS 
internal energy, and Cheat is the energy due to shock heat- 
ing Cheat = e — e cold- 



V. NUMERICAL METHODS 

Duez et al. and Shibata and Sekiguchi ^} have in- 
dependently developed new codes to evolve magnetized 
fluids in dynamical spacetimes by solving the Einstein- 
Maxwell-MHD system of equations self-consistcntly. 
Both codes evolve the Einstein field equations without 
approximation, and both use high-resolution shock cap- 
turing techniques to track the MHD fluid. Several tests 
have been performed with these codes, including MHD 
shocks, nonlinear MHD wave propagation, magnetized 
Bondi accretion, MHD waves induced by linear gravi- 
tational waves, and magnetized accretion onto a neu- 
tron star. Details of our techniques for evolving the 
Einstein-Maxwell-MHD system as well as tests can be 
found in |lfiL Il7|. In this paper, we have performed sev- 
eral simulations for identical initial data using both codes 
and found that the results are essentially the same. 

The simulations presented in this paper assume ax- 
ial and equatorial symmetry. We evolve only the x-z 
plane [a (2+l]_dimensional problem]. We adopt the car- 
toon method [3!j for evolving the BSSN equations, and 
use cylindrical coordinates for evolving the induction and 
MHD equations. In this scheme, the coordinate x is iden- 
tified with the cylindrical radius w, and the y-direction 
corresponds to the azimuthal direction. For example, for 
any vector V\ V x = V™ , and V y = zoV*. 

When black holes appear in our simulations, we avoid 
the singularity by using black hole excision. This tech- 
nique involves removing from the grid a region inside the 
event horizon which contains the spacetime singularity. 
Rather than evolving inside this region, boundary condi- 
tions are placed on the fields immediately outside. For 
details on our excision techniques, see |32l 1331 140| . 

As in many hydrodynamic simulations, we add a ten- 
uous, uniform-density "atmosphere" to cover the com- 
putational grid outside the star. For stars A, Bl, and 
B2, the rest-mass density in the atmosphere is set to 
p a = 10~ 7 /9 max (0), where p m ax(0) is the initial maximum 
rest-mass density. The initial pressure in the atmosphere 
is set to the cold polytropic value (P = Kp\\). If the 
density in a given grid cell drops below p a after an evo- 
lution step, we simply set p — p a . We also impose limits 
on the pressure in order to prevent negative values of the 
internal energy and to prevent spurious heating of the 
atmosphere. In particular, if the pressure drops below 
Praia = 0-5Kp r , we set P = P m i n ; similarly, if P rises 
above P max = 10Kp r , we set P = P max . Our main re- 
sults are not sensitive to the adopted (small) value of p a ; 
similarly for P min and P max . 

Due to the hybrid EOS, we found that a different atmo- 
sphere scheme is appropriate when evolving star C. For 
this case, we choose p a = 10 9 g/cm 3 ps lCT 6 /9 max (0). The 
specific internal energy e of the atmosphere is set to be 
ifi(100/9 a ) ri_1 /(ri — 1) = £ m in- If the value of e becomes 
smaller than this value, we artificially set e — e m in- We 
also limit the maximum value of e as 30e co id; if the value 
of e exceeds this value, we artificially set e = 30e co id- 
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VI. NUMERICAL RESULTS 
A. Star A 

We have performed simulations on star A with fixed 
initial field strength (C = 2.5 x 10~ 3 ). We use a uniform 
grid with size (N,N) in cylindrical coordinates (w, z), 
which covers the region [0, L) in each direction. We have 
performed simulations with L = 4P oq and 5P oq and 
found that the results depend only weakly on L. In the 
following, we present results with L = 4.5i? oq . For star A, 
P oq = 4.5M = 18.6 km(M/2.8M Q ). To check the con- 
vergence of our numerical results, we perform simulations 
with four different grid resolutions: N = 250, 300, 400 
and 500. Unless otherwise stated, all results presented 
in the following subsections are from the simulation data 
with resolution N — 500. We will first describe the gen- 
eral features of the evolution and then discuss the effects 
of resolution, the behavior of the various components of 
the energy, and the excision evolution. 

1. General features of the evolution 

Figure ^ shows the snapshots of density contours and 
poloidal magnetic field lines (lines of constant A v ) in the 
meridional plane. Figure|3]shows the snapshots of three- 
dimensional (3D) magnetic field lines projected onto the 
equatorial plane. 

In the early phase of the evolution, the frozen-in 
poloidal magnetic fields lines are wound up by the differ- 
entially rotating matter, creating a toroidal field which 
grows linearly in time (see Fig. [3 and Fig. HJI) with the 
growth rate predicted by Eq. (JJJ. When the magnetic 
field becomes sufficiently strong, magnetic stresses act 
back on the fluid, causing a redistribution of angular 
momentum. The core of the star contracts while the 
outer layers expand. At t > 6P C , the effect of the MRI 
becomes evident, as shown in Figs. 0J; and where 
we see that the maximum value of \B X \ (= \B m \) sud- 
denly increases, growing exponentially for a short period 
(about one e- folding). We find that the MRI first oc- 
curs in the outer layers of the star near the equatorial 
plane. This is consistent with the linear analysis, as 
Eq. (|12|l together with star A's angular velocity profile 
gives a shorter tyim near the outer part of the star. The 
effect of the MRI can be seen in Fig. ^ where we see 
that the poloidal field lines are distorted. The growth of 
the central density slows down once l-B^max and |P y | max 
(= |tuP v | max = |P T | max ) saturate at t ~ 20P C . This 
may be caused by MRI-induced turbulence redistribut- 
ing some of the angular momentum to slow down the 
contraction of the core. The amplitude of the toroidal 
field begins to decrease after t > 20P C ~ (see Figs. El 
and 0Ji) and the core of the star becomes less differen- 
tially rotating (Fig.EJ. This is consistent with the results 
of 21], which predict that the magnetic field growth by 
magnetic winding should saturate after an Alfven time, 



the magnetic energy having grown to an appreciable frac- 
tion of the initial rotational kinetic energy. 

The combined effects of magnetic braking and MRI 
eventually trigger gravitational collapse to a black hole at 
t ~ 66P C w 36(M/2.8Af Q ) ms when an apparent horizon 
forms. A collimated magnetic field forms near the polar 
region at this time (see Fig. [TJ. However, a substantial 
amount of toroidal field is still present (see Fig- ■ With- 
out black hole excision, the simulation becomes inaccu- 
rate soon after the formation of the apparent horizon be- 
cause of grid stretching. To follow the subsequent evolu- 
tion, a simple excision technique is employed 33, 40]. We 
are able to track the evolution for another 300M w 8P C . 
We find that not all the matter promptly falls into the 
black hole. The system settles down to a quasiequilib- 
rium state consisting of a black hole surrounded by a hot 
torus and a collimated magnetic field near the polar re- 
gion (see the panels corresponding to time t — 74.6P C 
in Fig. 2]). The irreducible mass of the black hole is 
about 0.9M and the rest-mass of the torus is about 0.1M 
(Fig.0). We estimate that J/M 2 ~ 0.8 for the final black 
hole. This system is a promising central engine for the 
short-hard gamma-ray bursts (see Sec. IVlEl and 20]). 



2. Resolution study 

Four simulations were performed with different resolu- 
tions (see Fig. N =250, 300, 400 and 500. We find 
that the results converge approximately when N > 400. 
On the other hand, results are far from convergent for 
N < 300. For example, |-B x | max is much smaller at lower 
resolutions than for runs with higher resolutions, and 
the growth rate of l-B^lmax is underestimated. Hence, 
the effect of MRI, which is responsible for the growth of 
I B x | max , is not computed accurately for low resolutions. 
This is because the wavelength of the fastest growing 
MRI mode is not well-resolved for low resolutions. We 
find that we need a resolution A/A max < 0.14 (N > 400) 
in order to resolve the MRI modes. The straight dashed 
line in Fig.^Ji corresponds to the linear growth rate pre- 
dicted by Eq. Q. This slope agrees with the actual 
growth of l-B^lmax in the early (magnetic winding) phase 
of the simulation, but as back-reaction (magnetic brak- 
ing) becomes important, the toroidal field begins to sat- 
urate. 

Figure [5] shows the onset of the MRI in more detail for 
the two highest resolutions. Also shown is an approxi- 
mate fit to the growth rate (the short dashed line) . This 
line shows that the perturbation grows approximately as 
SB X oc e wt , where u> w 0.18/P C . This is a somewhat 
lower rate than that predicted from linear theory, which 
gives w m ax ~ 1/Pc, where w max corresponds to the fastest 
growing MRI mode. This discrepancy could be due to 
the fact that the linear analysis is inaccurate by a signif- 
icant factor. One drawback of the linear analysis is the 
assumption of Newtonian gravity, but star A is highly 
relativistic. In addition, the linear analysis treats the 
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FIG. 1: Snapshots of rest-mass density contours and poloidal magnetic field lines for star A at selected times. The first and 
third rows show snapshots of the rest-mass density contours and velocity vectors on the meridional plane. The second and 
fourth rows show the corresponding field lines (lines of constant A v ) for the poloidal magnetic field at the same times. The 
density contours are drawn for po/p m ax(0) = io~°- 36l ~ - 09 (i — 0-10), where p max (0) is the maximum rest-mass density at t = 0. 
The field lines are drawn for A v = A^min + (A v<max — A Vim i n )i /20 (i = 1-19), where A VtYaax and ^li^min are the maximum and 
minimum values of A v , respectively, at the given time. The thick solid (red) curves denote the apparent horizon. In the last 
panel, the field lines are terminated inside the black hole at the excision boundary. 



MRI as a purely local phenomenon, assuming a uniform 
background state over length scales much longer than the 
wavelengths of the perturbations. However, since the ex- 
pected A max is only one order of magnitude smaller than 



the initial equatorial radius, these assumptions may lead 
to significant discrepancies between the predicted and ac- 
tual properties of the MRI. 
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FIG. 2: Snapshots of the projected 3D magnetic field lines for star A (see Appcndix[X]for details) at selected times. Only three 
lines are drawn in each panel to prevent overcrowding of field lines. 
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FIG. 3: Angular velocity profiles for star A at selected times 
(corresponding to the times in Fig. 0. The last two profiles 
correspond to the moment of excision and a late time in the 
excision run. 



3. Evolution of the energies vs. time 




t/P 

» c 

FIG. 4: Evolution of the central rest-mass density p c , central 
lapse ct c , and maximum values of \B X \ and \B V \. IB 1 ] max and 
|B H |max are plotted in units of y/ p m ax(0). The solid (red), 
long-dashed (green), dot-dash (blue), and dotted (black) 
curves denote the results with A?=250, 300, 400, and 500 re- 
spectively. The dashed (cyan) line in (d) represents the pre- 
dicted linear growth of l-B^lmax at early times from Eq. J7J. 



Figure[B]shows the evolution of various energies defined 
in Sec. IIV Dl We see that the magnetic energy E^m re- 
mains small throughout the entire evolution, even though 
the magnetic field drives the secular evolution. The grav- 
itational potential energy W and the adiabatic part of 
the internal energy -Ei n t, a d change the most, which re- 
sults from the drastic change in the configuration of the 
star. The rotational kinetic energy decreases substan- 
tially before the core collapses, presumably because the 
bulk of the mass of the star rotates slower than at t = 0. 
A substantial amount of heat (£"hoat) is also generated 
by shocks. 



4- Evolution with excision 

Soon after the formation of the apparent horizon, the 
simulation becomes inaccurate due to grid stretching and 
an excision technique is required to follow the subsequent 
evolution. During the excision evolution, we track the ir- 
reducible mass of the black hole by computing the area of 
the apparent horizon ^4ah and using M\„ w \J Aar/IQtt. 
The irreducible mass and the total rest mass outside the 
apparent horizon are shown in Fig. [7| The total ADM 
mass of the final state system, consisting of a BH sur- 
rounded by a massive accretion torus, is well defined. In 
contrast, there is no rigorous definition for the mass of 
the black hole itself. To obtain a rough estimate, we 
proceed as follows. First, the angular momentum of the 
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FIG. 5: Evolution of [-B^lmax plotted in the same units as in 
Fig.[I]for the two highest resolution runs of star A. The dot- 
dashed (blue), and dotted (black) curves denote the results 
with iV=400 and 500, respectively. The dashed (cyan) line 
represents an approximate slope w = 0.18/P C for the expo- 
nential growth rate of the MRI, 5B X oc e ut . 
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FIG. 6: Components of the energy vs. time for star A. All 
energies are normalized to the binding energy at t = 0, where 
the binding energy is defined as -Ebind = Mo — M. In the 
evolution, i?t>md should be nearly conserved. 
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FIG. 7: Evolution of the irreducible mass and the total rest- 
mass outside the apparent horizon. (Here, t-ah is the local 
coordinate radius of the apparent horizon.) 
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FIG. 8: L2 norms of the errors in the Hamiltonian (H.) and 
momentum constraints (A4 Z ) for star A. The long-dashed, ver- 
tical (magenta) line represents the initial time for the excision 
run. 

as in Eq. I|49|). Then to estimate the black hole mass, we 
use 

M hole » yjM*; + (J ho i c /2M irr ) 2 , (61) 

which is an approximate relation for the spacetime of 
our numerical simulation, but would be exact for a Kerr 
spacetime. We thus find ~ 0.9M, where M is the 

total ADM mass of the system, and Jhoie/-Mhoic ~ 0.8. 

The black hole grows at an initially rapid rate follow- 
ing its formation. However, the accretion rate Mq gradu- 
ally decreases and the black hole settles down to a quasi- 
equilibrium state. By the end of the simulation, Mq has 
decreased to a steady rate of w 0.01M /P C , giving an 
accretion timescale of ~ 10-20P C w 5-10 ms(M/2.8M ). 
Also, we find that the specific internal thermal energy in 
the torus near the surface is substantial because of shock 
heating. The possibility that this sort of system could 
give rise to a GRB is discussed in Section lYl El and [2fj| . 

5. Constraint violations 



black hole is computed from 

•/hole = J - Jmatter^ > r AH ) (59) 

where the angular momentum of the matter outside the 
horizon is given by 

Jmatter(r > r AH ) = I Sytfx , (60) 

JV,r>r A H 



We monitor the violation of Hamiltonian and momen- 
tum constraints during the evolution. Figure |S] shows 
the L2 norm of the constraints. We see that in the 
pre-excision phase, the violation of all constraints are a 
few x 10~ 3 . Prior to excision, the constraints are satisfied 
to better than 1%. This indicates that our numerical evo- 
lution data accurately satisfy the constraint equations. 
After excision, the constraint errors jump to ~ 10%, but 
they remain constant for > 300M w 8P C . We thus can 
track the evolution reliably for > 2800M in total, which 
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is a nontrivial feat for highly relativistic, nonvacuum, and 
dynamical spacetime simulations. 
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FIG. 9: Selected parameters plotted against scaled time 
(t/tA) for evolutions of star A with four different magnetic 
field strengths: C = 1.25 x 1CT 3 (solid red lines), C = 
2.5 x 10~ 3 (green long-dashed lines), C = 5.0 x 10~ 3 (blue 
short-dashed lines), and C = 10 -2 (black dotted lines). All 
runs were performed with the same resolution (400 2 zones 
with outer boundaries at 20M) . When plotted against scaled 
time, the curves line up at early times (t sJ 0.5tA = H-Pc) 
when the evolution is dominated by magnetic winding. 



B. Star A, comparison of different values of C 

In order to test the scaling of our results for different 
values of the initial magnetic field strength, we have ex- 
amined three other values of C in addition to the value of 
2.5 x 1CP 3 chosen for the results of Section lVI Al Namely, 
we consider C = {1.25,2.5,5.0,10} x 10" 3 , and the re- 
sults are shown in Figs. |!|] and For the portion of the 
simulations in which magnetic winding dominates, the 
behavior is expected to scale with the Alfven time [2l| . 
In other words, the same profiles should be seen for the 
same value of t/tA- (The Alven time is inversely pro- 
portional to the magnetic field strength and hence pro- 
portional to C -1 / 2 .) From Fig. it is evident that this 
scaling holds very well for the toroidal field and for the 
central density and lapse, while t < OTi^. After the 
toroidal field saturates, the evolution is driven mainly by 
the MRI, which does not scale with the Alfven time. The 
scaling also does not hold during the collapse phase, when 
the evolution is no longer quasi-stationary. Though the 
scaling breaks down at late times in these simulations, 
the qualitative outcome is the same in all cases. 

The behavior of l-B^lmax for these four different values 
of C is shown in Fig.^3 The sudden sharp rise of l-B^lmax 




FIG. 10: Maximum value of \B X \ plotted vs. t/P c for evolu- 
tions of star A with four different magnetic field strengths. 
The line styles correspond to the same values of C as in 
Fig-El The behavior of l-B^lmax is dominated by the effects of 
the MRI and thus does not scale with the Alfven time. The 
curves corresponding to the two highest values of C (dotted 
and dashed) terminate at the time when the star collapses. 



signals the onset of the MRI, and the approximate agree- 
ment of the slopes for different values of C indicates that 
the exponential growth rate of the MRI does not depend 
on the initial magnetic field strength (as expected from 
the linear analysis). In cases with a very weak initial 
magnetic field, turbulence induced by the MRI may be- 
come important much earlier than the effects of magnetic 
braking, since the timcscale for the growth of the MRI 
does not depend on the initial magnetic field strength. 
In this case, the scaling with tA would not hold during 
any phase of the evolution. However, since both the MRI 
and magnetic braking lead to similar angular momentum 
transfer, the qualitative outcome may again be the same. 



C. Star Bl 

Here, we present results for the evolution of star Bl 
with C = 2.5 x 10~ 3 . This run was performed with resolu- 
tion 400 2 and outer boundaries located at 4.5i? (36.4M). 
Since this model is not hypermassive, the redistribution 
of angular momentum through MHD effects will not lead 
to collapse. However, since this star is ultraspinning and 
angular momentum is conserved in axisymmetric space- 
times, it cannot relax to a uniform rotation state every- 
where unless a significant amount of angular momentum 
can be dumped to the magnetic field. We find that this 
model simply seeks out a magnetized equilibrium state 
which consists of a fairly uniformly rotating core sur- 
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FIG. 11: Evolution of central rest-mass density p c , central 
lapse a c , maximum values of \B X \ and \B V \ for star Bl. The 
magnetic fields l-B^lmax and |-B y | ma x are plotted in units of 
\/Vmax(0). Note that the lines become fairly horizontal at 
late times, indicating that an equilibrium has been reached. 
The dashed line in (c) represents an approximate slope of 
uj — (0.37/P c ) for the exponential growth rate of the MRI, 
8B X oc e . The dashed line in (d) represents the predicted 
linear growth of |-B y | m ax computed from Eq. 0. 



rounded by a differentially rotating torus. This is similar 
to the final state we found in [T2 for the same model 
when evolved with shear viscosity. 

Figure lTTl presents the evolution of some relevant quan- 
tities for this case. From the central density and lapse, it 
is evident that the star has settled into a more compact 
equilibrium configuration. This is consistent with the 
expectation that magnetic braking should transfer angu- 
lar momentum from the core to the outer layers. A brief 
episode of poloidal magnetic field growth due to the MRI 
is indicated by the plot of IP 31 1 max in Figure ITU The in- 
stability saturates and quickly dies away [4l| , leaving the 
strength of the poloidal field largely unchanged. Early in 
the evolution, the maximum value of the toroidal com- 
ponent \B V \ rises due to magnetic winding. This growth 
saturates at t ~ 10P C ~ 0.5£a . We note, however, that 
the toroidal magnetic field is non-zero in the final equi- 
librium state, though it is no longer growing due to mag- 
netic winding. The accuracy of the spacetime evolution is 
demonstrated by Fig. El which shows that the Hamilto- 
nian and momentum constraint errors remain very small 
throughout the simulation. 

Snapshots of the evolution in the x-z plane are shown 
in Fig. El The density contours for times t — through 
25.0P C show that angular momentum redistribution leads 
to the formation of a more compact star surrounded by 
a torus. At t — 10P C , the distortions of the magnetic 



FIG. 12: L2 norms of the errors in Hamiltonian and momen- 
tum constraints during the evolution of star Bl. 



field lines due to the MRI are clearly visible. As the 
disk expands, magnetic field lines attached to this low 
density material open outward, eventually leading to the 
field structure seen in the last four times shown in Fig.fTTjl 
for which some field lines are still confined inside the star 
while others have become somewhat collimated along the 
z-axis. For t J> 35P C , the density contours and poloidal 
magnetic field lines change very little, indicating that the 
system has reached an equilibrium state which is quite 
different from the initial state. The effects of magnetic 
braking in this case are demonstrated by the series of 
snapshots in Fig. El which is analogous to Fig. [5] The 
field lines become very tightly wound for t ~ 10P C and 
relax at later times. However, a significant toroidal field 
persists at late times when the system has essentially 
settled down to a final state. 

In order to understand the behavior of this case, we 
plot in Fig. El the degree of differential rotation AO, 
defined as follows: 

J (VI 2 ) - (O) 2 

AO = V{ [ - } , (62) 

where the angular brackets refer to density weighted av- 
erages ((/) = J d 3 xp*f /Mo) and (O) is the average an- 
gular velocity at t = 0. Rather than approaching zero 
at late times, this quantity approaches a roughly con- 
stant value. Thus, the equilibrium final state still has 
significant differential rotation. The evolution of the an- 
gular velocity profile for star Bl is shown in Fig. for 
the equatorial plane. Figures. El an d El suggest that 
the final state consists of a fairly uniformly rotating core 
surrounded by a differentially rotating torus. However, 
this differential rotation no longer winds up the magnetic 
field lines (i.e., the toroidal field strength does not grow). 
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FIG. 13: Snapshots of density contours and poloidal magnetic field lines for star Bl. The first and third rows show snapshots 
of the rest-mass density contours and velocity vectors on the meridional plane. The second and fourth rows show the field lines 
(lines of constant A v ) for the poloidal magnetic field at the same times as the first and third rows. The density contours are 
drawn for p /Pmax(0) = icr - 36 *" - 09 (i = 0-10). The field lines are drawn for A v = A VtJain + {A v , w - ^, roin )i/20 (i = 1-19), 
where yl^max and A^^n are the maximum and minimum values of A v respectively at the given time. Note that the field lines 
and the density contours show little change for t > 35P C , indicating that the star has settled down to an equilibrium state. 



This is because the rotation profile has adjusted so that 
is approximately constant along magnetic field lines. 
This is demonstrated in Fig. 1171 which shows that 



at late times. Since the rotation profile is adapted to 
the magnetic field structure, a stationary final state is 
reached which allows differential rotation and a nonzero 
toroidal field. 



(63) 



Since the final state is still differentially rotating and 
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FIG. 14: Snapshots of the projected 3D magnetic field lines for star Bl. 
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FIG. 15: Evolution of the degree of differential rotation AQ 
for star Bl. At late times, AQ approaches a non-zero constant 
value. This shows that the final equilibrium state of star Bl 
is still differentially rotating. 



FIG. 17: Evolution of {\B J djQ\) (normalized to unity at t = 
0). Note that this quantity drops toward zero at late time, 
indicating that the star is driven to a differentially rotating 
equilibrium state in which Q is constant along the magnetic 
field lines. 
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FIG. 16: Angular velocity profiles at selected times (corre- 
sponding to the times in Fig. 1131 for star Bl. 
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FIG. 18: Components of the energy vs. time for star Bl. All 
energies are normalized to the binding energy at t = 0. In the 
evolution, E hind should be nearly conserved. 



is threaded with magnetic fields, this configuration must 
be checked for the presence of the MRI. From the linear 
(and local) analysis discussed in Sectional we found that 
the predicted wavelength for the fastest growing mode 
is ~ 2 — 3M at late times, whereas the radius of the 
final star is ~ 6M. (Since the local analysis does not 
take into account gradients in the vertical direction, it 



is qualitative at best in this regime. However, this does 
suggest that A ma x ~ R-) Thus, the magnetic field is no 
longer weak, and the MRI is likely suppressed. This is 
corroborated by the fact that we do not see any rapid 
magnetic field growth at late times. 

Figure El shows the evolution of various energies. As 
in the case of star A, the magnetic energy -Eem shows a 
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FIG. 19: Evolution of magnetic energy Eem for star Bl. The 
energy is normalized by the initial rotational kinetic energy 
of the star, T rot (0). 



much smaller change in amplitude than T rot , Ei ntac i an d 
W. These results are very different from those found 
in 0, ^3 an d from those of the star B2 evolution (see 
the next subsection), where the change in Eem is com- 
parable to the change in T ro t- This is probably because 
star Bl is ultraspinning, in contrast to the "normal" 
models in 0, an d star B2. Here the seed magnetic 
field in star Bl causes a substantial change (on a secular 
Alfven timescale) in the structure of the star. The en- 
ergies -Emt.ad, ?r t and W readjust to the values of the 
new configuration, which is significantly different from 
the initial state. On the other hand, star B2 and the 
models studied in 0, 22 show little or no change in the 
density profile. As a result, a decrease in T rot results in 
an increase in Eem- 

Figure ^3 shows the evolution of the magnetic energy 
Eem, normalized to the initial rotational kinetic energy 
of the star, T rot (0). The value of EEM/T ro t(0) rises from 
its initial value of 6.7 x 1CP 4 to a peak of ~ 0.06 (the cor- 
responding field strength is about 90 times the initial field 
strength) mainly due to magnetic braking. Then it grad- 
ually decreases to the equilibrium value of 0.014. The 
final magnetic field strength |Bfl na i| is about 4.5 times 
the initial value. In cgs units, we find that for the initial 
field considered here, 

ISflnall ~ 10 17 {^f) G . (64) 

This field is comparable to the field strength of a magne- 
tar. Since the strength of the initial seed magnetic field 
is much smaller than the strength when it saturates, it 
is possible that the final equilibrium state will be the 
same even if the initial seed field is much smaller than 
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FIG. 20: Angular velocity profiles in the equatorial plane for 
star B2 at times t = [thick solid (black) line], t = 8.1P C ss 
It a [dashed (red) line], t = 32.5P C ~ 4£a [long dashed (green) 
line], and t = 46.5P C 5.8tA [dotted (blue) line]. At late time 
(t > 30P C ~ 4£a), the bulk of the star is nearly uniformly 
rotating. 



the present value. If this is true, a new-born neutron star 
with mass and angular momentum distribution similar to 
star Bl is likely to end up as a magnetar due to MHD 
processes. 



D. Star B2 

Both stars Bl and B2 are nonhypermassive. However, 
star Bl is ultraspinning, whereas B2 is normal. We evolve 
this star with a seed magnetic field strength C — 2.5 x 
10~ 3 . Our simulation shows that this star evolves to 
a uniformly rotating configuration with little structural 
change (see Figs. 1201 and l21f) . 

Figure l2Tl shows the density contours and poloidal mag- 
netic field at the initial time (t = 0) and at t = 46.5P C « 
5.8^- We see that the density profile of the star does not 
change appreciably. This is not surprising since the main 
effect of the MHD processes is to redistribute the angu- 
lar momentum inside the star. However, the rotational 
kinetic energy of star B2 is not very large (the initial 
T/|W| = 0.040). Hence, the change of the centrifugal 
force inside the star as a result of angular momentum 
transport does not disturb the initial equilibrium signif- 
icantly, unlike the cases of stars A, Bl and C (see the 
next section). 

Figure 1^21 shows the evolution of various energy com- 
ponents. Unlike stars A, Bl and C (see Fig. l28[) . the mag- 
netic energy Eem and rotational kinetic energy T rot show 
the largest fractional variations. The adiabatic internal 
energy Ei nt . a d and gravitational potential energy W have 
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FIG. 21: Snapshots of the rest-mass density contours and 
poloidal magnetic field lines for star B2 at times t = and 
t = 46.5P C - The first row shows snapshots of the rest-mass 
density contours on the meridional plane. The second row 
shows the corresponding field lines for the poloidal magnetic 
field at the same times. The density contours are drawn for 



Po/Pmax(0) = 10 



(i = 0-10), where p ma x(0) is the 



maximum rest-mass density at t = 0. The field lines are 



drawn for A^, — A v 



Oi/15 (i = 1-14), 



values of A v , respectively, at the given time. The meridional 
components of the velocity (which are zero initially) at t — 
46.5P C are very small and so are not shown here. 



very small fractional changes. This is because the con- 
figuration of the star does not deviate significantly from 
the initial equilibrium (see Fig. I21f> . A large fraction of 
the growth of magnetic energy comes from the rotational 
kinetic ener gy ( see Fig.l2r>[). This is similar to the results 
reported in pllli^ . 
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FIG. 22: Components of the energy vs. time for star B2. All 
energies are normalized to the binding energy at t — 0. Some 
quantities are normalized by an additional numerical factor 
(as indicated) to ease visualization. 
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FIG. 23: Change of T TOt and £7 mag vs. time. Here, ST IC 
Trot -T rot (0) and 8E mag — ^mag ^mag (0). 



We next demonstrate that the same qualitative fea- 
tures of the MHD-induced hypermassive collapse dis- 
cussed in Section IVl Al are also present with a more real- 
istic EOS. To do this, we evolve star C, which was con- 
structed using the hybrid EOS described in Section ITTT1 
The ADM mass of this star is 2.65M Q , which is 17% 
larger than the mass limit of a rigidly rotating neutron 
star for the adopted hybrid EOS. We choose an initial 
magnetic field with C = 7.1 x 10~ 3 as the fiducial model. 
In this case, the maximum strength of the magnetic field 
is ~ 5 x 10 16 G. The computational domain is [0, L] in 



the x- and z-directions, with L = 5R C 



54km. We 



performed the same evolution with resolutions N = 501, 
601, and 751 to check convergence. 

Snapshots of the evolution at eight selected times are 
shown in Fig. El Figure 123 shows the evolution of the 
maximum density, central lapse, maximum values of \B X \ 
and \B y \ for the three values of N, indicating approxi- 
mate convergence. The maximum values of \B X \ and \B y \ 
increase as the value of N is increased. This is a natural 
consequence of the fact that the profile of the magnetic 
field is better resolved with increasing N. 
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FIG. 24: The same as Fig.Qbut for star C. The contours for the first and third rows are drawn for po = lo 15-0 4 * g/cm 3 (i — 0- 
9). In the last two panels, curves with po = 10 11 g/cm 3 (solid curves) and with po = 5 x 10 10 g/cm 3 (dotted curves) are 
also drawn. The circular arc near the bottom-left corner in last three panels denotes an apparent horizon. The second 
and fourth rows show the poloidal magnetic field lines at the corresponding times. The solid contour curves are drawn for 
Ai/, — 0.8(1 — Q.li)Atp,max,o (i = 0-9) and the dotted curves are for A v — 0.08(1 — 0.2i)A v , m ax,o (i = 1-4). Here, A v = A v ,max,o 
is the maximum value of A v at t = 0. Note that the outer computational boundary in this simulation is located at x w 54 km 
and z fa 54 km, and that P c « 0.2 ms. The results with N = 601 are shown here. 



As in the case of star A, the early phase of the evolution 
(t 13P C ) is dominated by magnetic winding. The lin- 
ear growth then saturates and the subsequent evolution 
is dominated by the MRI (see the snapshots at t = 11.5P C 
in Fig. [2] in which a clear distortion of the poloidal mag- 
netic field lines is seen for 1 km <J x Si 4 km and z 3 
km) . Soon after the onset of the MRI, the outer layers of 



the stellar envelope are blown off (see the snapshots for 
t/P c = 11.5-31.1). This explosion causes an expansion 
and redistribution of the magnetic field lines. Eventu- 
ally, the removal of angular momentum from the central 
regions by the MRI results in collapse and black hole 
formation at t ~ 33P C . 

The winding up of the toroidal magnetic field leads to 
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FIG. 25: Evolution of the maximum rest-mass density p ma x, 
central lapse a c , and maximum values of \B X \ and \B V \ for 
star C. I-B^lmax and |-B a | max are plotted in units of \/ p max (0). 
The dashed (green), long-dashed (blue), and solid (black) 
curves denote the results with resolutions of N = 501, 601, 
and 751 respectively. The dotted lines in (c) and (d) corre- 
spond to an exponential growth rate of 1/3P C and predicted 
linear growth of l-B^max at early times from Eq. respec- 
tively. 



strong, inhomogeneous magnetic pressure. The toroidal 
field is primarily generated in regions where the initial 
poloidal magnetic field has a significant radial compo- 
nent. Thus, material at high latitudes gains a high 
magnetic pressure at early times in the evolution. In 
Fig. contour curves for the ratio of the magnetic pres- 
sure P mag = b 2 /2 to the gas pressure P are shown for 
t < 11.5P C . It is seen that the region around x ~ 5 km 
and z ~ 4 km has the maximum ratio P m &g/P, and the 
initial seed magnetic field is roughly radial in this region 
(see the first snapshots of Fig. 124(1 . This region of strong 
magnetic pressure beneath the surface of the star is sub- 
ject to the effects of magnetic buoyancy [43, H3, and 
toroidal magnetic field lines suddenly emerge from inside 
the HMNS, propelling material outward in the explosion 
(see the snapshots of Fig. 1261 for t/P c <; 8). This be- 
havior may be due to the interchange instability |4g ■ A 
similar magnetic buoyancy phenomenon is also observed 
in star A. However, unlike star C, the magnetic buoyancy 
does not cause an explosion in star A's outer layers. 

The time scale for the rearrangement of the field 
and the fluid due to buoyancy is approximately 
the same as that of the convection instability, and 
hence, of order r buoy - {R 2 H/GM) 1 / 2 {p/ ^p) 1 ' 2 or ~ 



(R 2 H/GM) 1 r 2 {P/AP mas ) 1 ' 2 [4jj where R and H are the 
equatorial radius and scale height of the inhomogeneity 
of magnetic pressure and Ap/p is the degree of inhomo- 



geneity of the density due to the inhomogeneity of the 
magnetic pressure AP mag . For the outer layers of star C, 
we find H ~ 2 km, and AP mag is approximately equal to 
the magnetic pressure P mag . Since (P/Pnag) 1 ^ 2 ~ c s /Va, 
we have 

rbuoy ~ {R 2 H/GM) 1 / 2 c s /V A ~t A c,/(GM/H)V 2 

~ 0At A ~ 0.9 ms , (65) 

which is comparable to the Alfven time scale. Indeed, 
this churning of field lines and fluid due to magnetic buoy- 
ancy seems to begin as soon as the toroidal field is wound 
up to a significant strength. 

The formation of the black hole is accompanied by the 
formation of a torus (see the last three snapshots of Fig. 
124(1 . To follow the growth of the black hole due to accre- 
tion, the subsequent evolution of the system is computed 
with excision. Since the torus is magnetized, turbulent 
motion is induced which transports angular momentum 
outward in the accretion torus and encourages the accre- 
tion of matter onto the black hole. 

In the accretion torus, the magnetic fields have a strong 
radial component (see, e.g., the snapshots at t = 44.5P C ). 
This is because, during the formation of the black hole, 
some material in the envelope of the HMNS is ejected in 
the radial direction (see the snapshot at t/P c = 31.1), 
enhancing the radial magnetic field. The ejected matter 
soon falls onto the accretion torus, and compresses the 
magnetic fields. This process leads to a strong magnetic 
field in the accretion torus. (The typical value of P mag /P 
is 10 3 -10 4 near the surface of the accretion torus.) As a 
result, material from high latitudes (which is originally 
blown away from the HMNS as a wind during the col- 
lapse) does not fall toward the equatorial plane, but col- 
lides with the surface of the torus, and then falls into 
the black hole along the surface of the torus (see the vec- 
tor fields at t = 44.5P C ). Hence, accretion occurs along 
high latitudes as well as along the equatorial plane. This 
scenario for black hole accretion is slightly different from 
those presented, e.g., in [5(1153]. We also note that the 
last three panels of Fig. 1241 show that the density of the 
accretion torus gradually decreases, indicating that ac- 
cretion is quite rapid in the first ~ 15P C after the forma- 
tion of the black hole. 

For t J> 45P C , the accretion relaxes to a steady rate 
M ~ 5 X W- 4 M /P C ~ 5M /s. The final state con- 
sists of a rotating black hole surrounded by a hot torus 
undergoing quasistationary accretion. At t = 50P C , the 
irreducible mass of the black hole is M; rr w 0.9M, while 
the torus consists of ~ 1% of the original rest mass and 
~ 4% of the original angular momentum of the system 
(see Fig. 127(1 . During the simulation, ~ 1% of the total 
rest mass and ~ 5% of the total angular momentum es- 
cape from the computational domain through outflows. 
Following the same calculations as in Sec. IVI A 41 we es- 
timate the mass and spin parameter of the black hole at 
t « 50P C to be M ho ie « 0.98M and J ho i e /M£ ole w 0.75. 

Figure |2H1 shows the evolution of the various energies 
defined in Sec. HVDl The magnetic energy Pem reaches 
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FIG. 26: Contour curves for the ratio of the magnetic pressure P mag to the gas pressure P for t < 11.5P C . The contour curves 
are drawn for P mag /P = (f , ma g /P)max X 10" 01 (magenta), 10" 5 (red), 10" 1 (blue), and lO" 05 ' (i = 3,4) (black). The time 
t is indicated for each snapshot. 



a value of at most 8% of the binding energy (£bmd) 
throughout the entire evolution, even though the mag- 
netic field drives the secular evolution. The gravitational 
potential energy W and the adiabatic part of the inter- 
nal energy -Ei n t,ad change the most, which results from 
the drastic contraction of the stellar core. The fraction 
of Sint.ad is 60-70% larger than that for star A. This is 
simply due to the fact that star C is more compact. The 
rotational kinetic energy is nearly constant. A substan- 
tial amount of heat (-Bheat) is generated by shocks. When 
the apparent horizon first appears, this heat is ~ 1% of 
the rest-mass energy (i.e., ~ 5 x 10 52 ergs). Most of the 
heat is swallowed by the black hole, but a substantial 
fraction remains in the accretion torus (see below). Fi- 
nally, the binding energy decreases by « 7% by the end 
of the pre-excision evolution. This is mainly due to the 
violation of approximate conservation of the ADM mass 
by w 0.5% and to the escape of ~ 1% of the mass from 
outer boundaries. 



The internal energy in the torus corresponds to a typi- 
cal thermal energy per nucleon of approximately 10 2 MeV 
[2fj|. giving an equivalent temperature T Ril-2 xlO 11 
K for the density ~ 10 10 -10 12 g/cm 3 if the assumed 
components are free nucleons, ultra-relativistic electrons, 
positrons, neutrinos, and thermal radiation 44] . The 
opacity to neutrinos inside the torus (considering only 
neutrino absorption and scattering interactions with nu- 
cleons) is Q 



K ~ 7X 10_17 cn^g- 1 . (66) 



Because of its high temperature and density, the torus 
is optically thick to neutrinos. Thus, the neutrino lumi- 
nosity is estimated 43] as L„ ~ nR 2 F, where R is the 
typical radius of the emission zone and F is the flux from 
the neutrinosphere. In the diffusion limit, F is approxi- 
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FIG. 27: Post-excision evolution of star C. The irreducible 
mass Mi rr of the black hole and the rest mass of the torus 
surrounding the black hole settle down to their quasiequilib- 
rium values at late times. The results with N = 601 are 
shown. 
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FIG. 28: Components of the energy vs. time for star C. All 
energies are normalized to the binding energy at t = 0, where 
the binding energy is defined as Ehmd = Mo — M. In the 
evolution, i?bmd should be nearly conserved. 
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where a is the Stefan-Boltzmann constant, N v is the 
number of thermal neutrino species, taken as 3, and £ is 
the surface density of the torus ~ 10 17 -10 18 g/cm 2 . We 
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FIG. 29: Evolution of the maximum rest-mass density p max , 
central lapse a c , and maximum values of |B H |/y p ma x(0) as 
a function of t/tA for star C with three values of C. The 
solid (blue), dashed (green), and long-dashed (red) curves 
correspond to results with C = 1.55 x 10 -2 , 7.1 x 10~ 3 , and 
3.8 x 10 -3 , respectively. The grid size is N = 601 for all 
cases. The dotted line in (c) corresponds to the predicted 
linear growth of |B H | m ax at early times from Eq. 



then obtain 



L v ~ 2 x 10 w ergs/s 



T 



R 



10 km 



10 11 KJ \10 17 g/cm 2 



(68) 



This luminosity will be present for the total duration of 
the accretion, ~ 10 ms. Since the torus has a geomet- 
rically thick structure, a substantial fraction of neutri- 
nos are emitted toward the rotation axis, leading to en- 
hanced neutrino-antineutrino pair annihilation along the 
axis. The pair annihilation could produce a relativistic 
fireball since the baryon density near the rotation axis 
is much lower than that in the torus. Furthermore, the 
luminosity is expected to have a strong time-variability 
because of the turbulent nature of the torus. Therefore, 
this massive and hot torus has many favorable properties 
which may explain a short GRB of energy ~ 10 48 -10 49 
ergsPy). This possibility was explored by Shibata et al. 

Two other simulations are performed (with N — 601) 
for different values of the initial magnetic field strength: 
C = 3.8 x 10~ 3 and 1.55 x 10~ 2 . In Fig. HO we show 
the evolution of the maximum density, central lapse, and 
maximum value of \B V \ as a function of t/tA- For star C, 
we find t A /P c = 8.94, 13.4, and 17.9 for C = 1.55 x 10~ 2 , 
7.1 x 10~ 3 , and 3.8 x 10~ 3 , respectively. Figure 1291 shows 
that the scaling relationship holds for t/tA ^ 1 as in 
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FIG. 30: Evolution of the maximum values of \B X \/^/ p ma x(0) 
as a function of t/P c for star C with three values of C. The 
solid (blue), dashed (green), and long-dashed (red) curves 
denote the results with C = 1.55 x 1CT 2 , 7.1 x 1CT 3 , and 
3.8 x 1CT 3 , respectively. The grid size is TV" = 601 for all 
cases. The dotted line segment corresponds to an exponen- 
tial growth rate of 1/3P C . 

Fig. El The scaling breaks down when t/tA > 1, indi- 
cating that the MRI and other effects such as magnetic 
buoyancy determine the evolution of the system. 

In Fig. 1301 we show the evolution of the maximum 
value of | B x | max as a function of t/P c for three values 
of C. As in Fig. QJJI the sudden exponential growth sig- 
nals the onset of the MRI, and the approximate agree- 
ment of the growth rate for different values of C indi- 
cates that the exponential growth rate of the MRI does 
not depend on the initial magnetic field strength. Af- 
ter the exponential growth, the magnitude of l-B^lmax 
remains roughly constant until the dynamical collapse 
occurs. During this phase before collapse, the angular 
momentum is transported outward gradually by the tur- 
bulence. The duration for this angular momentum trans- 
port is ~ 15P C irrespective of the value of C as long as 
C <; 3.8 x 10~ 3 . This indicates that the angular momen- 
tum is transported by a mechanism independent of the 
initial magnetic field strength (probably the turbulent 
transport associated with the MRI). 

We note that the collapse time of the HMNS reported 
here depends slightly on the parameters of the atmo- 
sphere, although the timescales for growth of the mag- 
netic field due to winding and the MRI do not. This 
is inevitable since, just before the collapse, the HMNS 
is only marginally stable against a quasiradial instabil- 
ity, and thus, a slight increase in the atmospheric mass- 
energy sensitively shortens the collapse time. 

We have also studied models with masses slightly dif- 
ferent from that of star C presented here. We find that 
the mass of the resulting torus varies significantly. For 



more massive stars, the torus mass is smaller. This is 
probably due to the fact that the star collapses sooner 
and hence there is less time for outward angular momen- 
tum transport. For a sufficiently large mass, the resulting 
torus mass is smaller than 0.1% of the total mass, which 
is probably too small for the system to trigger a short 
GRB. On the other hand, less massive stellar models re- 
sult in larger torus masses. This result is interesting since 
it might explain the variety of short GRBs. The details 
of this study will be reported in a future paper. 



VII. DISCUSSION AND CONCLUSIONS 

We have discussed in detail the evolution of magne- 
tized HMNSs as first reported in [l^. |2C| . In addition, 
we have performed simulations of two differentially ro- 
tating, but nonhypermassive, neutron star models with 
the same initial magnetic field geometry. These simula- 
tions have revealed a rich variety of behavior with pos- 
sible implications for astrophysically interesting systems 
such as binary neutron star remnants, nascent neutron 
stars, and GRBs, where magnetic fields and strong grav- 
ity both play important roles. 

The two hypermassive models considered in this study, 
stars A and C, both collapse to BHs due to the influence 
of the initially poloidal, seed magnetic field. The early 
phase of evolution for both models is dominated by mag- 
netic winding. As the strength of the toroidal magnetic 
field grows, the resulting magnetic stress begins to trans- 
port angular momentum from rapidly moving fluid ele- 
ments in the inner region to the more slowly moving fluid 
elements in the outer layers. During this magnetic brak- 
ing phase, the inner regions of the stars undergo quasis- 
tationary contraction, while the outer layers expand and 
begin to form a low-density torus. 

The winding of the magnetic field proceeds until the 
back-reaction on the fluid becomes strong enough that 
the growth of the toroidal field ceases. This happens 
after roughly one Alfven time. After several rotation pe- 
riods, we also see the effects of the MRI. Plots of the 
poloidal magnetic field lines display perturbations with 
wavelengths similar to A max (the wavelength of the fastest 
growing mode estimated from the linear analysis). These 
perturbations first appear in the outer layers of the star, 
which is consistent with the linear analysis. In order to 
diagnose the sudden local growth of the poloidal mag- 
netic field due to the MRI, we track the maximum value 
of \B X \ on the grid. We found that IP^Imax grows expo- 
nentially at a rate which does not depend on the strength 
of the initial magnetic field, in accord with the proper- 
ties of the MRI. However, the growth rate observed in our 
numerical simulations differs significantly from that pre- 
dicted by the linear analysis. This is due probably to the 
fact that (a) the linear analysis is based on Newtonian 
gravity, but the models we study here are highly relativis- 
tic, and/or (b) the small MRI wavelength assumption in 
the analysis might not be applicable to our magnetic field 
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configuration. 

The nonlinear outcome of the MRI is turbulence, 
and this turbulence leads to further angular momentum 
transport. Eventually, the inner cores of stars A and C 
become unstable and collapse to BHs. Surrounding the 
BHs, significant amounts of material remain in magne- 
tized tori which have been heated considerably by shocks 
resulting from the turbulent motions of the fluid. This 
final state consisting of a BH surrounded by a massive, 
hot accretion disk may be capable of producing highly 
relativistic outflows and a fireball (either through v — v 
annihilation or MHD processes) and is hence a promis- 
ing candidate for the central engine of short-hard GRBs. 
This model predicts that such GRBs should accompany 
a burst of gravitational radiation and neutrino emission 
from the HMNS delayed collapse. 

The behavior of the nonhypermassive, ultraspinning 
star Bl under the influence of a seed magnetic field is 
quite different. Magnetic braking and the MRI operate in 
this model as well, leading to a mild contraction of the in- 
ner core and the expansion of the outer layers into a high 
angular momentum torus-like structure. The final state 
consists of a fairly uniformly rotating core surrounded 
by a differentially rotating torus. The remaining differ- 
ential rotation does not shear the magnetic field lines (i.e. 
(|i? J <9jT2|) approaches zero in the final state), so that the 
toroidal field settles down. We find that this configura- 
tion is not subject to the MRI, probably because it is sup- 
pressed by the strong magnetic field (the MRI wavelength 
is comparable to the size of the star, and the standard 
local linear analysis breaks down in this regime) . The ro- 
tation state of the final configuration naturally depends 
on the geometry of the initial magnetic field. On the 
other hand, the normal star B2 simply evolves to a uni- 
formly rotating configuration. 

Two issues in particular warrant further study. The 
first is the scaling behavior of our solutions. We begin 
our simulations with a seed magnetic field which, though 
far too weak to be dynamically important, may be sig- 
nificantly larger than magnetic fields present in HMNSs 
formed through stellar collapse or a binary neutron star 
merger. We have demonstrated that, by varying the 
strength of the initial magnetic field through a factor of 
~ 3 (See Fig. El, our evolution obeys the expected scaling 
during the magnetic winding phase, and the qualitative 
outcome of the simulations remains the same. However, 
since the MRI grows on a timescale ~ few x P c regard- 
less of the initial magnetic field strength, it is possible 
that, for very weak initial fields, the effects of the MRI 
could dominate the evolution long before the effects of 
magnetic braking become important. In this case, the 
scaling of our numerical results with the Alfven time (rel- 
evant for magnetic winding) may break down. The rel- 
ative importance of magnetic winding and the MRI for 
different seed field strengths deserves further study. Un- 
fortunately, the wavelength of the fastest growing MRI 
mode becomes very difficult to resolve numerically as the 
strength of the initial magnetic field decreases. However, 



our results seem to indicate that magnetic braking and 
MRI-induced turbulence have similar effects in magne- 
tized HMNSs. Thus, the qualitative features of the evo- 
lutions described here may also be present for HMNSs 
with much weaker initial seed fields. 

Another issue which warrants further study concerns 
the effects on our evolutions of relaxing the axisymme- 
try assumption. Rapidly and differentially rotating neu- 
tron stars may be subject to bar and/or one-armed spi- 
ral mode instabilities which could affect the dynamics 
(though star A was shown in 0, 0] to be stable against 
such instabilities, at least on dynamical timescales). Ad- 
ditionally, the development of the MRI in 2D differs from 
the 3D case Turbulence arises and persists more 

readily in 3D due to the lack of symmetry. More specif- 
ically, according to the axisymmetric anti-dynamo theo- 
rem , sustained growth of the magnetic field energy 
is not possible through axisymmetric turbulence. This 
phenomenon has been demonstrated by numerical simu- 
lations 47] . However, McKinney and Gammie [5(| have 
performed axisymmetric simulations of magnetized tori 
accreting onto Kerr BHs and have found good quan- 
titative agreement with the 3D results of De Villiers 
and Hawley f° r the global quantities E/M Q and 
J/Mq which are the rates of total energy and an- 
gular momentum falling into the horizon, normalized by 
the accretion rate. Though simulations in full 3D will 
eventually be necessary to capture the full behavior of 
magnetized HMNSs, the 2D results presented here likely 
provide (at least) a good qualitative picture. 
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APPENDIX A: DRAWING MAGNETIC FIELD 
LINES 



The vector potential A 1 is related to the magnetic field 
B i by B l = n^ k djA k , where e a ^ 5 is the Levi-Civita 
tensor. It is easy to show that in axisymmetry, the 
poloidal components of a magnetic field (£? ro and B z ) 
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are determined by A v alone as follows: 
vj^jB™ = ~d z A v , 



(Al) 
(A2) 



Poloidal magnetic field lines are two-dimensional curves 
on which dzu/dz — B^ / B z — ~d z A v /d^A v . Hence we 
have 

dA v = (d^A v )dw + (d z A v )dz = (A3) 

on the curves. This means that contours of constant A v 
are the poloidal magnetic field lines. All the poloidal field 
lines shown in this paper are drawn by the contours of 
Am. There are two ways of calculating Am. One method 
is to integrate Eqs. l|Aljl and l|A2(l . The other method is 
to evolve A v according to the equation 



dtAp = w^j{v z B™ - v™B z 



(A4) 



This equation can be derived from Eqs. (35) and (42) 
of |H3. 

In order to show the toroidal component of the mag- 
netic field, we draw field lines projected onto the equato- 
rial plane. To do this, we first choose three points 'cun) 
(j =1, 2, 3) in the equatorial plane so that at the given 
time, 



when there is no apparent horizon in the time slice. If 
there is an apparent horizon, we set z m j n = if ra/fj) > 

0.5r AH and z min = ^J0.25r 2 AH - w 2 {j) if tu (i) < 0.5r AH - 

Here tah is the coordinate radius of the apparent hori- 
zon. Next we integrate the equations 



dx 



d\ 
dX 



dz 



U) 



dX 



= B x (x ij) ,y {j) ,z U) ) , 
= B v {x {j) ,y {j) ,z {j) ) , 
= B z {x ( j ) ,y {j) ,z l j ) ) , 



(A6) 
(A7) 
(A8) 



with the initial locations (xj,yj, Zj) given by: 

"2(j-l)7rl 



X(J)\\=0 

V(j)\x=o 

ZQj)\\=0 



G7(j) COS 



^(j) sin 

^min 



2(j - IK 



(.7 = 1,2,3) 



(A9) 



Am(tU(^), z m j n ) — Am, m i n -t- (A(p iQiax Am ;m i n )j/4 , (A5) 

where A VtVnax and Am !m i n are the maximum and mini- 

mum values of A v at the given time [H3 |. and z m [ n = (y/ x 2 + y 2 , z) 



Here A serves as a parameter of the 3D curves. The 
integration is terminated when the curve goes beyond 
the boundary of the grid. The projected field lines are 
the trajectories (x(j) (A), (A)) traced out by A. On 
the other hand, the poloidal field lines determined by 
the contours of A v are equivalent to the trajectories of 
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